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1. Introduction 

Consider the classical problem of computing the expected value of a real function 
/ of the d-variate random variable X as a linear combination of its values f{z) at a 
finite set points z S I? C M''. The general quadrature problem is: determine classes 
of functions f : M.'^ ^ M., finite set of n nodes I? C K'' and positive weights {Xz}z^v 
such that 

(1) E(/(X)) = / fix) d\{x) = /(^)^- 

where A is the probability distribution of the random variable X. In the simplest 
univariate case, d = 1, the set D is the set of zeros of a node polynomial, e.g. the 
n-th orthogonal polynomial for A, see e.g. [H Sec. 1.4]. Not much is known in the 
multivariate case, unless the set of nodes is a product of one-dimensional set. 

A similar setting appears in statistical Design of Experiment (DoE) where one 
considers a finite set of treatments V and experimental outputs as funtion of the 
treatment. The set of treatments and the set of nodes are both described efficiently 
as zeros of systems of polinomial equations, i.e. as what is called in Commutative 
Algebra a 0- dimensional variety. Such a framework has been used in a systematic 
way in the literature on Algebraic Statistics using the tools of modern Computa- 
tional Commutative Algebra, see e.g. [8] and [5]. In such studies the set V is called 
a design and the ajfine structure of the ring of real functions on T) is analyzed in 
detail because it represents the set of real responses to treatments in V. However, 
the euclidean structure, such as the computation of mean values, is missing in the 
algebraic setting. In algebric design of experiment the computation of mean val- 
ues has been obtained by considering very special sets called factorial designs, e.g. 
{+1, —1}'', see e.g. [3]. Note that {-1-1, —1} is the zero set of the polynomial — 1. 

The purpose of the present paper is to discuss how both worlds can be treated 
together by considering orthogonal polynomials. In particular, we consider algo- 
rithms from the world of Commutative Algebra for the cubature problem in ([T|) by 
mixing tools from elementary orthogonal polynomial theory and from probability. 
Viceversa, the formula ([Ij provides an interesting interpretation of the equation in 
the RHS as expected value. 

We proceed in sequence by increased degree of difficulty. 

In Section [5] we consider the univariate case and take A to admit an orthogonal 
system of polynomials. Let g{x) — Jldepl^ ~ ^^"^ univariate division given a 
polynomial p there exist unique q and r such that p — q g+r and r has degree smaller 
than the number of points in V. Furthermore r can be written as X^deP '''{d)ld{x) 
where Id is the Lagrange polynomial for d Cz V. Then we show that 

(1) the expected values of p and r coincide if and only if the n-coefficients of 
the Fourier expansion of q with respect to the orthogonal polynomials is 
zero 

(2) the weights are the expected values of the Lagrange polynomials. 
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When A is a standard Gaussian probability law and V the zero set of the n-th 
Hermite polynomial, the application of Stein-Markov theory premits a representa- 
tion of some Hermite polynomials, including those of degree 2n — 1, as sum of an 
element in the polynomial ideal generated by the roots of the Hermite polynomial 
of degree n and of a reminder, suggests a folding of multivariate polynomials over 
a finite set of points. 

The point is to describe a ring structure of the space generated by Hermite 
polynomials up to a certain order because it is essentially the aliasing on functions 
induced by limiting observations to T). The particular structure of the recurrence 
relationship for Hermite polynomials makes this possibile and we suspect that the 
study of the ring structure over 2? for other systems of orthogonal polynomials will 
require different tools from those we use here. 

This result implies a system of equations in Theorem 13.51 (extented to the mul- 
tidimensional case in Section [SJ which gives an implicitly description of design and 
weights via two polynomial equations. We envisage applicability of this in the choice 
of V for suitable classes of functions but have not developed this here. The case 
when the design is a proper subset of the zero set of the n-th Hermite polynomial 
is developed in Section H) 

Section [6] contains our most general set-up: we restrict to ourselves to product 
probability measures on but consider any set of n distinct points in M''. Then 
a Buchberger-Moller type of algorithm is provided that works exclusively with or- 
thogonal polynomials. It gives a generating set of the vanishing ideal of T) written 
in terms of orthogonal polynomials. Then it is applied to compute the coefficients 
of ((H in the Fourier expansion of the interpolatory polynomial at T) of the function 
whose expectation is wanted. Of course it will is of interest to determine generali- 
sations of our results to the cases where A is not a product measure and still admits 
an orthogonal system of polynomials. 

1.1. Basic commutative algebra. We start with some notation on polynomials: 
'R.[x\ is the ring of polynomials with real coefficients and in the d- variables (or 
indeterminate) x = (xi, . . . , Xd)] for a — (ai, . . . , Ud) a d-dimensional vector with 
non-negative integer entries, indicates a monomial; -<r indicates 

a term-ordering on the monomials in M.[x] i.e., a total well ordering on Z>q such 
that x" -<r x^ if divides x^ . If d = 1 there is only one term ordering of the 
monomials, which is definitely not the case for d > 2. Design of product form 
share some features with the dimension one. Because of that term orders are not 
much used in standard quadrature theory. We will see that refining the division 
partial order to a proper term-order is actually relevant in some cases, but not in 
all multivariate cases. 

The total degree of the monomials x" is X)iLi '^i- The symbol M.[x\k indicates the 
set of polynomials of total degree k and R[x]<fc the vector space of all polynomials of 
at most total degree k. Let P be a finite set of distinct points in K'', A a probability 
measure over and X a real- valued random variable with probability distribution 
A so that the expected value of the random variable f{X) is E(/(X)) = / f{x)d\[x). 

Given a term ordering -<t, for p e R[a;] there exists unique g and r e M[x] such 
that 



r is a linear combination of less monomials than points in P, and r is smaller than 
g with respect to -<r- It is a consequence of the multivariate polynomial division, 
see e.g. [5J Sec. 2.3]. Furthermore r{x) can be written uniquely as 



(2) 



p{x) — g{x) + r{x), g{x) — if x ^ T), 



(3) 
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where Iz is the indicator polynomial of the point z in V, i.e for a; G 2? it is lz{x) = 1 
\i X = z and lz{x) ={) \i x ^ z. Equation Q follows from the fact that {^^ : z G 2?} 
is a M-vector space basis of C{T>), the set of real valued functions over V. 

The expected value of the random polynomial function p{X) with respect of A 

is 

E {p{X)) - E {g{X)) + E {r{X)) - E {g{X)) + ^ p{z) E 

zev 

by linearity. In this paper we discuss classes of polynomials p and design points T) 
for which 

E{p{x)) = Y,pwmx)) 

zev 

equivalent ly E (.g(X)) = 0. 

The polynomials g and r in (0) are fundamental in the applications of algebraic 
geometry to finite spaces. In multivariate polynomial division if /i, . . . , /f G M.[x\ 
form a Grohner basis with respect to -<r (see [H Cap. 2]) and generate the ideal of 
polynomial functions vanishing over 2?, then 

t 

g{x) = hi{x)fi{x), hi{x) G M[a;] not unique 

i=l 

and the unique r has its largest term in -<-r not divisible by the largest term of 
gi, i = I,. . . ,t. Moreover, monomials not divisible by the largest terms of fi, i — 
1, . . . ,t, form vector basis of monomial functions for the vector space £(2?) of real 
functions on 2?. Various general purpose softwares, including Maple, Mathematica, 
Matlab and computer algebra softwares, like CbCcA, Macaulay, Singular, allows 
manipulation with polynomial ideals, in particular can compute the reminder and 
the monomial basis. 

For all p G R[a;], the polynomial r above is referred to as reminder or normal 
form. It is often indicated with the symbol NFt-(p, {/i, . . . , /t}), or the shorter ver- 
sion NF(p), while (/i, . . . , ft) indicates the polynomial ideal generated by /i, . . . , /t. 

In one dimension, a Grobner basis reduces to a polynomial / vanishing over V 
and of degree n = |2?| and r satisfies three main properties: 

(1) r is a polynomial of degree less or equal to n — 1, 

(2) p{x) — g{x)+r{x) — q{x)f{x) + r{x) for a suitable q G M.[x\ and g, / G (/), 
and 

(3) r{x) — p{x) if X is such that f{x) — 0. 

In the next Section we are going to study the algebra of orthogonal polynomials in 
one variable. 

2. Orthogonal polynomials and their algebra 

In this section let d = 1 and 2? be the set of zeros of a polynomial which is 
orthogonal to the constant functions with respect to A. We recall the basics on 
orthogonal polynomials we use next, see e.g. '4_. 

Let / be a finite or infinite interval of R and A a positive measure over / such 
that all moments /ij = x^ d\{x), j = 0, 1, . . ., exist finite. In particular, each 
polynomial function is square integrable on / and the L^{X) a scalar product is 
defined the ring M.[x] by 

{f{x),gix))^= J^fix)gix)d\ix) 

We consider only A whose related inner product is definite positive, i.e. ||/|| > 
if / 7^ 0. In this case there is a unique infinite sequence of monic orthogonal 
polynomials with respect to A and we denote them ttq, tti, . . .. Furthermore we have 
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TTfc G K[a;]fc, ttq, . . . jTTfc form a real vector space basis of R[x]<k, T^k is orthogonal 
to all polynomials of total degree smaller than k and for p € M\x\ and n E Z>o 
there exists unique c„{p) G R, called n-th Fourier coefficient of p, such that p(a;) = 
J2n=o Cn(p)7rn(a;) and only a finite number of c„(p) are not zero. 
Since the inner product satisfies the shift property 

{xp{x),q{x))^ = {pix),xq{x))^, 

then the corresponding orthogonal polynomial system satisfies a three-term recur- 
rence relationship. More precisely, all orthogonal polynomial systems on the real 
line satisfy a three-term recurrence relationships. Conversely, Favard's theorem 
holds [9]. 

Theorem 2.1 (Favard's theorem). Lei7„,a„,/3„ be sequences of real numbers and 
for n > let 7r„+i(a;) = {jnX — a„)7r„(a;) — /3„7r„_i(a::) be defined recurrently with 
TTo{x) — 1, 7r_i(a;) = 0. The 7r„(a;), n = 0,1,... form a system of orthogonal 
polynomials if and only if 7„ ^ 0, q„ ^ and an"/nln-i > for all n > 0. // 
7„ = 1 for all n then the system is of monic orthogonal polynomials. 

In the momc case, ak — —, r- and = hold true and the 

(7rfe,7rfe) (7rfc_i, 7rfe„i) 

norm of 7r„ is computed from the /3's as ||7r„|P = /3„/3n-i ■ • - /Jo- For orthonormal 

polynomials TTfe = 7rfc/||7rfc|| the Christofjel-Darboux formulae hold 

E" , s~ /— 7f„(x)7r„_i(i) - 7f„_i(x)7r„(t) 
nk[X)TTk[t) = VPn : 
X — t 

(4) 



n-1 
fc=0 



Non-example 2.2. Inner products of the Sobolev type, namely {u,v)s — {u,v)\g + 
{u' , w') Ai ■ ■ ■ + ^^'^•')a, where \i are positive measures possibly having different 
support, do not satisfy the shift condition. Neither do the complex Hermitian inner 
products. 

Theorem 2.3. Let I? = {a; C M : 7r„(x) = 0} be the zero set of the n-th orthogonal 
polynomial with respect to the probability measure A. If p{x) = q(x)7r„(a;) +r(x), 
then 

lK{p{X)) = ^p(z)A^, if and only if Cn{q) =0. 

Remark 2.4. This theorem is a version of a well known result, see e.g. PI Sec. 1.4]. 

We include the proof to underline a particular form of the error in the quadrature 
formula, to be used again in the next Theorem and in Section 6. 

Proof. The zero set of 7r„, T> = {x : 7r„(x) = 0}, contains n distinct points. For 
a univariate polynomial p, we can write uniquely p{x) = q{x)'Kn{x) + r{x) with 
deg(r) < n and deg(g) = max{deg(p) — n,0}. Furthermore, the indicator func- 
tions in the expression r{x) — Lagrange polynomials for T>: 

-.— w X — UJ 

namely Izix) — Ylwev-w^z z E T>. Hence we have 



E (piX)) = E {q{X)iT^{X)) + p{z) E {U{X)) 

zev 

= c„((?)||7r„||| + 'Yp{z)\^. 



zev 

□ 
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A particular case of Theorem 12 .31 occurs if p has degree less than 2n. In this case 
q has degree at most n — 1 and c„(q) = 0. This shows that the quadrature rule with 
n nodes given by the zeros of 7r„ and weights {A^lzgu is a Gaussian quadrature 
rule and it is exact for all polynomial functions of degree smaller or equal to 2n—l. 
For notes on quadrature rules see for example jH Chapter 1]. 

Example 2.5 (Identification). For / polynomial of degree iV < 2n — 1 we can write 
/(^) — X]fc=o Cfc(/)7rfc(a;). The constant term is given by 

coif) = mix)) J2 /(^)^- 

zev 

and for all i such that N + i < 2n — 1 

Mldf) - EifiX)TT,iXj) = J2 fi^)^^{^)K 

In particular, if deg f = n—1 then all coefficients in the Fourier expansion of / can 
be computed exactly. 

In general for a polynomial of degree N possibly larger than 2n — 1, Theorem l2.3l 
gives the Fourier expansion of its reminder by 7r„ , indeed 

/(^)^«(^)A. = J2 NF(/7r,)(z)A. = E(NF(/(X)7r,(X))) = |k,||2c,(NF(/)) 

Theorem 12.61 below generalises Theorem 12.31 to a generic finite set of n dis- 
tinct points in R, say V. As above, the indicator function of z g I? is Izix) = 

-—^ X — w 

Yiw&v-w^z • di^) ~ Yizevi^ ~ ^) unique monic polynomial 

vanishing over T) and of degree n. Write a polynomial p S M.[x\ uniquely as 
pix) — q{x)g{x) + r(x) and consider the Fourier expansions of q and g: qix) = 
Y.t=Q<^kiq)T^kix) andg(x) =YTk=oCk{g)'^k{x). 

Theorem 2.6. With the above notation, E(p(X)) = '^z^DPiz)\z if and only if 
E^J^Cfc(g)cfc(5)||7rfc||2 =0. 

Proof. From 

+00 n 

p{x) = qix)gix) + r{x) = ^ Ckiq)Tikix) ^ Cj(g)7rj(x) + ^ piz)U 

k=0 j=0 z<£V 

we have 

+C30 n 

EipiX)) = ^^Cfc(g)c,(.9)E(^fe(XV,(X)) + 5]p(z)A. 

n 

= Yc.kiq)ckig)\\TTk\\i + Yp{z)\:, 
k=o zev 
and this proves the theorem. □ 

The condition we find in Theorem 12.61 is linear in the Fourier coefficients of q, 
which is found easily from / by polynomial division. The first Fourier coefficients 
of q appearing in the conditions of the theorem are determined by solving the system 
of linear equations 

(5) M [ck{q)]k=o,...,\v\-i = b(^)]fc=o....,|i3|-i 

where M = [Kk{z)]^^xi k=o ^^"^ design/evaluation matrix for the first jX'l 

orthogonal polynomials. 

Theorem 12.61 can be used in two ways at least. If p is known, the condition in 
the theorem can be checked to verify if the expected value of p can be determined 
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by Gaussian quadrature rule with nodes V and weights E ^ntoeP to5<tz ^ ^^-"^ 

z &T>. The Fourier coefficients of g can be computed analogously to those of q 
adapting Equation If p is unknown and p{x) = X^aPa^;" for a finite number of 
non-zero, unknown real numbers pa, Theorem 12.61 characterizes all the polynomials 
for which the Gaussian quadrature rule is exact, namely E (p{X)) = X^zeP 
Furthermore, the characterization is a linear expression in the unknown p^, where 
P(^) — X^Q^'c^"- This is because in Equation ([5]) the q{z) are linear combinations 
of the coefficients of p. 

In Section [3] we shall specialise our study to Hermite polynomials, while in Sec- 
tion 16.11 we shall generalise Theorem 12.61 to higher dimension. To conclude this 
section, we observe that the remainder r admits an interpretation as a projection. 

Proposition 2.7. Let p{x) € W^x] and write p{x) — q{x)'Kn{x) + r{x) where r has 
degree less than n. Then q is the unique polynomial such that p — qiTn is orthogonal 
to all TTfn with m > n. 

Proof. As r = p — qiTn has degree at most n — 1, it can be written as a linear 
combination of ttq, tti, . . . , 7r„_i and is orthogonal to tt^ for all m > n. Let qi 
and q2 G M.[x] such that both p — qiHn and p — q2i^n are orthogonal to tt^ for all 
m > n. Now (qi — q2)'^n has degree not smaller than n, hence it is orthogonal 
to ttq, . . . , 7r„_i. Necessarily it is qi — q2 — 0, because it is orthogonal to all tt^, 
fc = 0, .... □ 

Example 2.8. Substituting the Fourier expansions of q and p in the orthogonality 
relation, we find that the m-th coefficient in the Fourier expansion of p can be 
written as 

E {p{X)Trrn{X)) - E ((Z(X)^„(X)^™(X)) 

^Cfc(p)E(7rfc(X)7r,„(X)) = ^c,(g)E(7r,(X)7r„(X)7r,„(X)) 



c™(p) I |7r„ 1 12 = ^ c, (g) E (tt, (X)7r„(X)7r„ (X)) 



For Hermite polynomials it can be simplified by e.g. using the product formula in 
Theorem 13.11 of Section [31 

3. Hermite polynomials 

There is another way to look at the algebra of orthogonal polynomials that we 
discuss here in the case of Hermite polynomials. The reference measure A is the 

normal distribution, dX{x) — w{x) dx, with w{x) — — , x £ R. 

V 27r 

3.1. Stein-Markov operators for standard normal distribution. For a real 
valued, differentiable function /, define 

Sf{x) = xf{x) - /(x) = -e-V2i. , 

d" = -; — , and consider Z ^ \. The following identity holds 
ax" 

(6) E (0(Z) (5'V(^)) = E {dJ'4){Z) tp{Z)) 

if 0, are such that linix^±oc, 4>{x)'4'{x)e^^ ^"^ — and are square integrable, see 
[21 Ch. V Lemma 1.3.2 and Proposition 2.2.3]). Polynomials satisfy these conditions 
and S is also called the Stein-Markov operator for the standard normal distribution. 
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Ho = 


1 


Hi{x) 


= X 


H2{x) 


— X 


Hsix) 


= x 


Hi{x) 


— X 


Hsix) 


= X 



The n-th Hcrmite polynomial can be defined as Hn{x) = (5"1. Direct computa- 
tion using d proves the following well-known facts 

(1) the first Hermite polynomials are 



- 3a; 
"4 -6x^ + 3 
lOx^ + 15x 

(2) Hnix) = (-l)"e^'/2^"(e-"^'/2) (Rodrigues' formula) 

n\ 

(3) d6—Sd ~ id from which the relationships dHn = ni7„_i, d"^Hn ~ — -Hn-m 

ml 

for m < n and the three-term recurrence relationship 
(7) Hn+i = xHn - nHn-i 

are deduced. 

(4) Hermite polynomials are orthogonal with respect to dX(x) = — dx. 

V 27r 

Indeed from Equation ([5]) we have E {Hn{Z)H„i{Z)) = nlSn^m where Sn,m = 

if n 7^ m and Sn,m — I ii n ^ m. 
Already we mentioned that {Hn{x) : n < d} spans M[a;]<ci and that Hn is 
orthogonal to any polynomial of degree different from n. The ring structure of the 
space generated by the Hermite polynomials is described in Theorem 13. II 



Theorem 3.1. The Fourier expansion of the product of H^Hn is 

M fk\ ., 

HkHn — Hn+k + [ i j \ { I *'^"+'i'-2i 



Proof. Note that {(j), = M{(j){Z)tp{Z)) is a scalar product on the obvious space 
and let rt < fc with Z ^ Af{0, 1) and ip, (j) square integrable functions for which 
identity ^ holds. Then 



' n 



= (i7„+fc,1/^)+^ ( . jfc(fc-l)...(fc-Z+l)(i/„+fc_2.,V) 

i=l 



□ 



Example 3.2 (Aliasing). As an application of Theorem l3.1l observe that the three- 
term recurrence relation for Hermite polynomials Equation ([7]) 



Hn+l — xHn — nHn-l 

evaluated on the zeros of Hn{x), say I?„, becomes Hn+i{x) = —nHn-i{x) where 

^n-l ,n+k TT 



indicates that equality holds for x £ 2?„. In general let Hn+k = X]?=o h^j^'^Hj 
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be the Fourier expansion of the normal form of Hn+k at I?„, where we simphfied 
the notation for the Fourier coefficients. Substitution in the product formula in 



lower order Hermite polynomials: 



Theorem 13.11 gives the formula to write ft," ' in terms of Fourier coefficients of 



Equating coefficients gives a closed formula 



nAk / \ /, 

v-^ I n\ I k 



In Table [T] the normal form of Hk+n with respect to _ff„ is written in terms of 
Hermite polynomials of degree smaller than n. For example, Hn+7i{x) — —n{n — 
l)(n — 2)Hn~3{x) + 3n_ff„_i(x) for those values of x such that Hn{x) — 0. 



k 


Hn+k = 


1 


-nHn-1 




2 


-n{n - l)Hn- 


-2 


3 


—n{n — l){n - 


^2)Hn-3+SnHn-l 


4 


—n(n — l)(n - 


^ 2)(n - 3)H„„4 + 8n{n - 1)F„_2 


5 


n! u 

(n-5)!-""-5 


f bnHn-i + 15n(n - l)(n - 2)i7„_3 


6 


n! u 

(n-6)!-""-6 


f 24n(n - l)(n - 2)(n - 3)i/„_4 + 10?i(n - l)(2n - b)Hn-2 



Table 1. Aliasing of Hn+k, fc = 1, . . . , 6 over V = {iJ„(a;) — 0} 



Example 3.3. Observe that if / has degree n + 1 equivalently fc = 1 then 

n-l 

/ = E ("'~)^'' + ^"(Z)-^" + Cn+l{f )Hn+l 

i=0 

n-2 

= E ^i(f)^i + - nCn+l{f)) Hn-1 

and all coefficients up to degree n — 2 are "clean" . 

We give another proof of Theorem 12.31 for Hermite polynomials. 

Corollary 3.4. Let Vn = {x : Hn{x) = 0} and p e R[x]. Let p{x) = q{x)Hn{x) + 
r(x) with, the degree of r smaller than n and let Z ^ A/'(0, 1). Then 

E {p{Z)) = ^ p{z)K if and only if E = 

zeVn 

-.— w X — zu 

with = E {lz{Z)) and U{x) = IltoePito^^z z-w ' ^ ^ 
Proof. From Equation ^ we have 

E(g(Z)i/„(Z)) =E(g(Z)(5"l) =E(dXZ)) 
Now by the same steps followed in the proof of Theorem 12.31 we conclude that 

E(p(Z))=E(d>(Z))+ p{z)K 

□ 
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3.2. Algebraic characterisation of the weights. Theorem 13.51 gives two poly- 
nomial equations whose zeros are the design points and the weights. The proof is 
based on the Christoffel-Darboux formulae, Equations (|4]). 

Theorem 3.5. Let P„ = {x : Hn{x) = 0}. 

(1) There exists only one polynomial A of degree n — 1 such that X{xk) — Xk 
for all k — 1, . . . ,n, 

(2) furthermore Xk = ^"^ n'^' ^n-ii^k)- Equivalently 

(3) the polynomial X satisfies 

Proof. The univariate polynomial A is the interpolatory polynomial of the values 
Afe's at the n distinct points in I?„ and hence it exists, unique of degree n — 1. To 
prove item 2., observe that for Hermite polynomials a„ — 0, /3„ = n, Hn{x) = 
Hn{x)/Vni and H'^ix) = ^/nHn-i{x)- Substitution in the Christoffel-Darboux 
formulae and evaluation at 2?„ = {xi, . . . , a;„} give 



(8) Hki^'^Hki^j) = if z ^ J ^'^(^O' - nHn-i{. 

k=0 k=0 

In matrix form Equations ([5]) becomes 

In^n dia.g{Hn-i{xi)'^ : i = 1, . . . ,n) 



Xif 



where H„ is the square matrix H„ = Hj{xi) and diag indicates 

L J z— l,,,,,n;j— 0,...,n — 1 

a diagonal matrix. Observe that ]HI„ is invertible and 

M-i = min'^ Amg{H-1^{x,) : z = 1, . . . ,n) 

Now, let / be a polynomial of degree at most n— 1, a typical remainder by division 
for iJ„, then f{x) = YJ]=o CjHj{x). Write / = ll„c where / = [/(a:j)]i=i,...,n and 
c — [cj]j. Furthermore note that 

c = H- 1/ - H^n-i dmg{H-^,{x,) : * = 1, . . . , n)/ 

= Win-' di&g{H~l^{xi)f{xi) : z = 1, . . . ,n) 



1 " 

(9) = - E H,{x,)f{^^)H-\{x.) 



n 

Apply this to the fc-th Lagrange polynomial, f{x) = lk{x), whose Fourier expansion 
is f{x) = J2]=o '^kjHj(x). Using lk{xi) = Sik in Equation obtain 

(10) - ^H,{xk)H-^,{xk) 

The expected value of lk{Z) is 

n-1 

Xk = E {lk{Z)) = Y Ck3 E (Hj{x)) = Cfco 
i=0 

Substitution in Equation ([TU| for j — Q gives Xk = ^H~'^^{xk) = ^"~"^^' iJ,7-i(^fc)- 
This holds for all fc = 1, . . . , n. 

Item 3. is a rewriting of the previous parts of the theorem because the first 
equation _ff„ {x) — states that only values of a; G are to be considered and the 
second equation is what we have just proven. □ 
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Item 2 in the theorem states that the weights are strictly positive. Furthermore, 
Theorem 12.31 apphcd to the constant polynomial p{x) = 1 shows that they sum to 
one. In other words, the mapping that associates z to Xz, z G "Dn, is a discrete 
probability density. Then Theorem 12.31 states that the expected value of the poly- 
nomial functions oi Z ^ A/'(0, 1) for which c„(g) = 0, is equal to the expected 
value of a discrete random variables X given by {X ~ Xk) = E(/fe(Z)) ~ A^, 
fc = 1, . . . , n 

n 

E(p(^))=^p(xfc)Afc = E„(p(^)) 
fc=i 

Example 3.6. For n = 3 the polynomial A in Theorem 13. 51 can be determined by- 
hand. For larger values of n an algorithm is provided in Section r3.3l The polynomial 
system to be considered is 

= H:i{x) = - 3a; 

2/3 = \{x)Hl = [Oq + Oix + e2X^){x^ - 1)2 

where A(a;) — 6o + Oix + 92X^. The degree of A(x)iJ| is reduced to 2 by using 
x^ — 3x 

(11) 2/3 = X{x)Hl =00 + OiAx + (90 + 4^2)^^ 

Coefficients in Equation PT|) are equated to give X{x) = ^ — 

In some situations, e.g. the design of an experimental plan or of a Gaussian 
quadrature rule, the exact computation of the weights might not be necessary 
and X{x) is all we need. When the explicit values of the weights are required, 
the computation has to be done outside a symbolic computation setting as we 
need to solve Hj,{x) = to get = {— \/3, 0, -\/3} and evaluate X{x) to find 
A_^ = A(-V3) = i = and Ao = A(0) - f . 

3.3. A code for the weighing polynomial. The polynomial A(a;) in Theorem l3.5l 
is called the weighing polynomial. Table 13.31 gives a code written in the specialised 
software for symbolic computation called CoCoA [1]. to compute the Fourier ex- 
pansion of A(a;) exploiting Theorem l3.5l 

Line 1 specifies the number of nodes N . Line 2 establishes that the working 
environment is a polynomial ring whose variables are the first (TV — 1)-Hermite 
polynomials plus an extra variable w which encodes the weighing polynomial; here 
it is convenient to work with a elimination term-ordering of Elim(w), so that 
the variable w will appear as least as possible. Lines 3, 4, 5 construct Hermite 
polynomials up-to-order N by using the recurrence relationships ([7]). Specifically 
they provide the expansion of Hj over with k < j for fc = 0,...,A^— 1. Line 6 
states than Hjy = HiHn^i — {N — 1)Hn-2 = 0, 'giving' the nodes of the quad- 
rature. Line 7 is the polynomial in the second equation in the system in Item 3 
of Theorem 13.51 and 'gives' the weights. In total there are N equations which are 
collected in an algebraic structure called an ideal whose Grobner bases [2] is com- 
puted in Line 8. In our application it is interesting that the Grobner bases contains 
a polynomial in which w appears alone as a term of degree one. Explicitly w in such 
polynomial provides the weighing polynomial written in terms of the first — 1 
Hermite polynomials. 

Line 9 in Table [331 gives the polynomial obtained for iV = 4, namely X(x) = 

^^j^ — ^-j§-- The nodes are ±y3±V6 and the values of the weights are 
showing that both nodes and weights are algebraic numbers but not rational num- 
bers. On a Mac OS X with an Intel Core 2 Duo processor (at 2.4 GHz) using 
CoCoA (release 4.7) the result is obtained for = 10 in Cpu time = 0.08, User 
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Line 1 N:=4; 

Line 2 Use R: :=Q [w,h[l . . (N-1)]] , Elim(w) ; 

Line 3 Eqs : = [h [2] -h [1] *h [1] +1] ; 

Line 4 For I: =3 To N-1 Do 

Line 5 Append(Eqs ,h [I] -h [1] *h [I-l] +(I-1) *h[I-2] ) EndFor; 

Line 6 Append(Eqs,h[l] *h[N-l] -(N-l)*h[N-2] ) ; 

Line 7 Append (Eqs, N*w*h [N-1] '2-Fact (N-1)) ; 

Lines J : =Ideal (Eqs) ; GB_J : =GBasis ( J) ; Last(GB_J); 

Line 9 3w + l/41i[2] - 5/4 

Table 2. Computation of the Fourier expansion of the weighing 
polynomial using Theorem 13.51 



time = 0; for TV = 20 in Cpu time = 38.40, User time = 38; for iV = 25 in Cpu 
time = 141.28, User time = 142 and for iV = 30 in Cpu time = 5132.71, User time 
= 5186 and gives a weighing polynomial of 22.349 characters. Observe that this 
computations can be done once for all and the results stored. 

4. Fractional design 

In this section we return to the case of general orthogonal polynomials, {7r„}„, 
and positive measure, dX. We assume that the nodes are a proper subset T of 
T^n — {x '■ T^nix) = 0} with a number of points m, < m < n. We work within two 
different settings, in one the ambient design I?„ is considered while in the other one 
it is not. 

Consider the indicator function of J- as subset of I?„, namely lj^{x) — I ii x (z T 
and if X £ I?„ \ J-. It can be represented by a polynomial of degree n because it 
is a function defined over I?„ [HIS] . Let p be a polynomial of degree at most n — 1 
so that the product p{x)ljr(x) is a polynomial of degree at most 2n — 1. Then from 
Theorem 12.31 we have 



E((pV)(X)) = Pi^)^z = E„ b(F)l^(y)) = E„ ip{Y)\Y € T) P„(y € F) 

where X is a random variable with probability law A and y is a discrete random 
variable taking value z E F with probability Vn{Y = z) = X^. The first equality 
follows from the fact that f{x)ljr{x) is zero for a; e 2? \ and the last equality 
from the definition of conditional expectation. 

Another approach is to consider the polynomial whose zeros are the elements of 
!F , say LiJj^{x) = W {x — z). Now consider the Lagrange polynomials for IF , namely 

lf{x) = n for z & F. 

w^z^ — w 

Lemma 4.1. Let J- C Vn- The Lagrange polynomial for z Cz IF is the remainder 
of the Lagrange polynomial for z € I?„ with respect to UJjr{x), namely 

l^Ax)^^^{lz{x),{u^{x))) 

Proof. There exists unique NF(Z2)(a;), polynomial of degree small than m, such 
that 

l,{x) = q{x)ujr{x) + ^¥{U){x) 

Furthermore, for a G we have lz{a) = NF(Z2)(a) — S^.a — The two 

polynomials l^ix) and NF(Z2)(x) have degree smaller than m and coincide on m 
points, by interpolation they must be equal. □ 



12 



CLAUDIA PASSING, GIOVANNI PISTGNE, AND EVA RICCGMAGNG 



For a polynomial p of degree N, write p{x) — q{x)ujjr(x) + r{x) with f(z) — r{z) 

if z e J" and r{x) — J2zeJ^Pi^)^'zi^)- Let q{x) — J2j=o"^ ^i'^ji^) ^^'^ ^t{x) = 
^™QCi7ri(x) as ujjr has degree m. Then 

(N-m m \ 

= ^OColkollA + ^lCl||7ri||;^ + . . . + 6(JV-m)AmC(JV —m)Am I \ '^{N—m)Am \ \ A 

where Af = E (NF(;2(X), {ujjr{X))), z £ T. 

Note that the error of the Gaussian quadrature rule, &oCo||7ro||A+&iCi||7ri||^ + . . .+ 
b(N-m)AmC(N-m)Am\\''^{N-m)/\m\\'\, IS linear in the Fourier coefficients bj, and also 
in the Fourier coefficients Cj relative to the node polynomial. This is generalised 
in Section 16.11 If the fraction J- coincides with the ambient design I?„ and hence 
contains n points and if p is a polynomial of degree at most 2n — 1 , then we obtain 
the well known result of zero error because (N ~ n) A n < n ~ 1 and the only 
non-zero Fourier coefficient of the node polynomial 7r„ is of order n. In general one 
should try to determine pairs of J- and sets of polynomials for which the absolute 
value of the errors is minimal. 



5. Higher dimension: zero set of orthogonal polynomials as design 

SUPPORT 

In this section we return to the higher dimensional set-up of Section 11.11 but 
we restrict ourselves to consider the product measure A'' = x^LiA and Xi, . . . ,Xd 
independent random variables each one of which is distributed according to the 
probability law A. As design we take a product grid of zeros of orthogonal poly- 
nomials with respect to A, more precisely our design points or interpolation nodes 
are 

= {x = {^ii . . . ,a;d) e M"* : TTm{xi) = TTnAx2) = ■■■ ^ T^n^i^d) = 0} 

where 7r„^ is the orthogonal polynomial with respect to A of degree n^. 

The Lagrange polynomial of the point y — (j/i, . . . ,yd) € I'ni,...,nd is defined as 
ly{xi, . . . ,Xd) = rifc^i ^yfc (^*:)' apex "*= indicates that ly^{xk) is the univariate 
Lagrange polynomial for yk € {xk : 7r„^ (xk) = 0} = Vn^ C M. 

The Span {ly '■ y € T^m,...,na)i is equal to the linear space generated by the mono- 
mials whose exponents lie on the integer grid {0, . . . , rii — 1} x . . . x {0, . . . , — 1}. 
Any polynomial / S can be written as 

d 

/(xi, ...,Xd)^'^qk{xi,.. ■ ,Xd)T^nk{xk) + r{xi,. ..,Xd) 

k=l 

where r is unique, its degree in the variable Xk is smaller than rifc, for A: = 1, . . . , d 
and belongs to that Span. 

The coefficients of the Fourier expansion of qk with respect to the variable Xk 
are functions of xi, . . . ,Xd but not of Xk- Let x-k denote the (d — l)-dimensional 
vector obtained from [xi, . . . ,Xd) removing the fc-th component and write 

d / +00 \ 

f{xi,...,Xd) = E ^Cj{qk)ix-k)Trjixk) TTn^ixk) + r{xi,...,Xd) 

k=l \j=0 J 

Only a finite number of Cj{qk){x-k) are not zero. 
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From the independence of Xi , . . . , Xn , the expected value of the Lagrange poly- 
nomial ly is 

d d 

{iy{x^, . . .,xa)) = n {iTM^)) = n 

k=l k=l 

where AJ!*" = E (lyl^{Xk)) is the expected value of a univariate random Lagrange 
polynomial as in the previous sections. 

Theorem 5.1. It holds 
E^. (/(Xi,...,Xrf)) = 

d 

fc=i {xi,...,x„)ev„-^,,,„^ 

Proof. The proof is very similar to that of Theorem l2.3l and we do it for d ~ 2 only. 
In a simpler notation the design is the nx m grid given by I?„m = {{x, y) '■ 7r„(a;) = 
= TTmiu)} and X and Y are independent random variables distributed according 
to A. The polynomial / is decomposed as 

f{x,y) = 

qi{x,y)nr,{x)+q2{x,y)TTn{y)+ ^ f{a,b)r^{x)lj;^{y)^ 

(a,6)eX'„,„ 

-t-oo -t-oo 

^Cj{qi){y)7rj{x) 7r„(a;) + ^ ((72)(a;)7rj- (y) 7r„(y)+ ^ f{a,b) l'2{x)l]!,\y) 

0=0 ]=0 (a,6)eD„.„ 

Taking expectation, using independence of X and Y and orthogonality of the tTj, 
we have 

{f{X,Y)) = 

E;,(c„((Zi)(y))|K„||2 +EA(c™(q2)W)||7r™||i+ J2 IMKK 



(a,f))eX'„,„ 



□ 



Note in the proof above that a sufficient condition for Ea {cn{qi)iY)) being zero 
is that / has degree in x smaller then 2n — 1, similarly for E^ {cmiq2){X)). We 
retrieve the well-known results that if for each i the degree in Xi of / is smaller 
than 2ni — 1, then 

E;,.(/(Xi,...,x,))= f{xi,...,xd)\:i...x:: 

(a:i,...,a:„)eX'„j...„^ 

In the Gaussian set-up, by Theorem 13.51 applied to each variable, weights and 
nodes satisfy the polynomial system 

Ai(a;i)i7„,_i(xi)2 







ni 



HnA^d) 

\d{xd)Hni^l{XdY 



= 



{nd - 1)! 



nd 



For the grid set-up of this section and for the Gaussian case, in analogy to 
Example 13.21 some Fourier coefficients of polynomials of low enough degree can be 
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determined exactly from the values of the polynomials on the grid points as shown 
in Example 15.21 below . 

Example 5.2. Consider a square grid of size n, and a polynomial / of degrees 
in X and in y smaller than n, the Hermite polynomials and the standard normal 
distribution. Then we can write 

ri-l 

As both the degree in x of fHk and the degree in y of fHh are smaller than 2n — 1, 
we have 

E{f{Zi,Z2)Hk{Xi)Hh{X2)) = Chk\\Hk{Xi)\\^\\Hh{X2)\\'' 

:T^i Y1 f{x,y)Hk{x)Hh{y)XxXy 

Note if / is the indicator function of a fraction T C then 

Cfe/i = -^jj^ ^ ^k{x)Hh{y)K\ with < /i, fc < n 

Example 15.31 deals with a general design and introduces the more general theory 
of Section [6] 

Example 5.3. Let F be the zero set of 

H2{x)^H2{y)=Q 
Hz{y)=Q 

Hi{x) (7J2(y)-2i/o) -0 

namely F is given by the five points (0,0), (±\/3, ±^3)- Write a polynomial 
/ £ K.[a;, y] s& f — ^ qiQi + r where r{x, y) = f{x, y) for (a;, y) £ F and r belongs 
to S^a.TL{Ho,Hi{x),Hi{y),Hi{x)Hi{y),H2[y)) = Span (l, x, y, xy, y^) . If, further- 
more, / is such that 

qi{x,y) = ao + aiHi{x) + a2Hi{y) + azHi{x)Hi{y) 
92 = 0i{x) + 92{x)H^{y) + e3ix)H2iy) 
q3 = a4 + ar^Hi{y) 

with a^Oj e M for i = 0, . . . , 5 and j = 1, . . . , 3, then E(.gi(Zi, ^2)9^(^1, Z2)) = 
for i = 1,2,3 and for Zi and Z2 independent normally distributed random variables. 
Write r as a linear combination of the indicator functions of the points in i.e. 
r{,x,y) = Y.(a.b)eJ' f(°'' ^)^{a,b)eJ^(^' v)- ^^'^^ indicator function l(a,b)f^T belongs to 
Spa,n{Ho, Hi{x), Hi{y), Hi{x)Hi{y), H2{y)) and they are 

2 1 , ^ 

l(o,o)e.F(a;,2/) = ^Hq - -^H2[y) 
hv-3,V-3)eA^^y) = + ^V3H,{x) + ^VSH^iy) + ^H,{x)H,{y) + ^H2{y) 

\Vz.-V-3)eA^^v) = Y^^o - + ^V3i/i(y) - l^H,{x)H,{y) + ^H2{y) 

\-V-z.V-^)eA^^y) = + ^^^i(^) - ^^iH,{y) - l^H,{x)H,{y) + ^^2(2/) 

l(-v^,-v^)e^(^'y) = T^^o - ^^H,{x) - ^V3H,{y) + l-H,{x)H,{y) + ^H2{y) 
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Their expected values are given by the i7o-coefficients. Furthermore, by linear- 
ity E(/(Zi,Z2)) -E(r(Zi,Z2)) =E(a,h)G^ IE(l(a,fa)6^(^i, ^2)) and we can 
conclude 

E(/(Zi,Z2)) =lE(r(Zi,Z2)) = 

^ /(0,0) , /(V3,V3) + /(73,-V3) + /(-V3,v^) + /(-73,-73) 
3 + 12 
The key points in Example 15.31 are 

(1) determine the class of polynomial functions for which E{gi{Zi, Z2)qi(Zi, Z2)) = 
and 

(2) determine the ifo-coeflicients of the indicator functions of the points in J- . 
In Section |6] we give algorithms to do this for any fraction J^. 

6. Higher dimension: general design support 

In the previous sections we considered particular designs whose sample points 
were zeros of orthogonal polynomials. In the Gaussian case we exploited the ring 
structure of the set of functions defined over the design in order to obtain recur- 
rence formula and to write Fourier coefhcients of higher order Hermite polynomials 
in terms of those of lower order Hermite polynomials (Example 13.21) . Also we 
deduced a system of polynomial equations whose solution consists the weights of 
a quadrature formula. The mathematical tool that allowed this is Equation ^ 
and the particular structure it implies for Hermite polynomials on the recurrence 
relation for general, orthogonal polynomials 

(12) T^k+iix) = (7fca; - ak)TTk{x) ~ l3kTTk-i{x) a; £ R 

with 7fc,afc 7^ and au^klk-i > 0. 

In this section we switch focus and consider a generic set of points in Mf^ as 
a design, or nodes for a cubature, and a generic set of orthogonal polynomials. 
We gain something and loose something. The essential computations are linear: 
such is the computation of a Grobner basis for a finite set of distinct points [7] ; the 
Buchberger Moller type of algorithm in Table[3]is based on finding solutions of linear 
systems of equations; in Section [5TT] we give a characterisation of polynomials with 
the same expected values which is a linear expression of some Fourier coefficients 
and a square free polynomial of degree two in a larger set of Fourier coefficients 
(see Equation [TO]) 

Given a set of points and a term-ordering the algorithm in Table [3] returns 
the reduced Grobner basis of the design ideal expressed as linear combination of 
orthogonal polynomial of low enough order. It does so directly; that is it computes 
the Grobner basis by working only in the space of orthogonal polynomials. 

We loose the equivalent of Theorem 13.11 for Hermite polynomials, in particular 
we do not know yet how to impose a ring structure on Span(7ro, . . . , 7r„) for generic 
orthogonal polynomials tt and we miss a general formula to write the product 7rfc7r„ 
as linear combination of tt^ with i — 0, . . . ,n A k^n + k, which is fundamental for 
the aliasing structure discussed for Hermite polynomials. 

For multivariate cubature formulae we refer e.g. to [10] which, together with [7], 
are basic references for this section. We are writing up in another manuscript our 
results on the degree of the cubature formula we obtain. For clarity we repeat 
some basics and notation. Let A be a one-dimensional probability measure and 
{7r„}„GZ>o be its associated orthogonal polynomial system. To a multi-index a = 
(ai, . . . , ad) & Z>Q associate the monomial x"^ ■ ■ ■ x^'', in short a;", and the product 
T^aiixi) ■ ■ ■ T^aaixd), in short TTaix). Notc that {TTajct are a system of orthogonal 
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polynomials for the product measure A"*. Theorem 16.11 describes the one-to-one 
correspondence between the a:"s and the 7rQ(a;)s. 

Theorem 6.1. (1) For d = 1 and k G Z>o, in the notation of Equation 
we have that 

k 

where cq{x^) = 1, c_i(.t") = ci(a;") — 0, and, for A; = 1, 2, ... , 
c-iix'') = Cfc+i(x'^') = 

7i-i 7i Ij+i 



1 



70 • • •7fe-i 

(2) For d > I, the monomial is a linear combination of tt^, with (3 < a 
component wise, and vice versa. In formulae 

(13) tTq = apx^ and x" = bjSTrp 

I3<a l3<a 

where (3 < a holds component wise. 

Proof. The proof of Item [T] is by induction and that of Item [5] follows by rearranging 
the coefficients in the product. See Appendix [71 □ 



Example 6.2. For ttj the j-th Hermite polynomial Hj, Item 1 of Theorem 16.11 
gives the well known result 

Cj{x'')^0 if A; + j is odd 

Cj{x'') = ^ ^ ^ (fc - j - 1)!! if k + j is even 

Direct application of Theorem 16.11 is cumbersome and we need only to charac- 
terise the polynomial functions for which the cubature formula is exact. So we 
proceed by another way. The finite set of distinct points 2? C M'* is associated to 
its vanishing polynomial ideal 

I{V) = {/ e R[x] : f{z) = for aU z e X>} 

Let LTa-{f) or LT{f) denote the largest term in a polynomial / with respect to 
a term-ordering a. Let [f(z)]zi£T> be the evaluation vector of the polynomial / at 
"D and for a finite set of polynomials G C M.[x] let [g{z)]z£V,g£G be the evaluation 
matrix whose columns are the evaluation vectors at I? of the polynomials in G. 

As mentioned at the end of Section [TTTl the space of real valued functions 

defined over 2? is a linear space and particularly important vector space bases can be 
constructed as follows. Let LT{I{V)) ^ (LT^lf) : f G If G is the cr-reduced 

Grobner basis of I{V), then LT{I{V)) = {LT„{f) : f &G). Now we can define two 
interesting vector space bases of C{V). Let L = {a € Z>q : x" ^ LT[I(V))}, then 
we define 

S = {x" : a e L} and OB = {tt^ : a C L} 

The sets L, B and OB depend on a. It is well known that if t G B and r divides t, 
then r € B; it follows that if a G i and (3 < a component wise then also /3 belongs 
to L and Tip to OB. For example for d = 2, let L = {(0, 0), (1,0), (0, 1), (2, 0)}, then 
B = {l,a;,2/,x2} and OB = {7ro(x)7ro(y), 7ri(a;)7ro(y), 7ro(x)7ri(y), 7r2(a::)7ro(?/)} = 
{1, 7ri(a;), TTi (y), 772(0:)}. Note that a induces a total ordering also on the orthog- 
onal polynomials: tTq tt^ if and only if a;" x^; analogously, with abuse of 
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notation, a <cr /? if and only if a;" for each a, /3 e ^>o- Further, given a < (3 

componentwise, since divides x^ and since x^~°' 1, we have <„ x^ , that 
is a <a P- 

Given a term-ordering cr, any g £ G can be uniquely written as its leading term, 
— LT(g), and tail which is a linear combination of terms in B preceding LT{g) 
in a, that is g = a:" + X];3gL i3<aa o-p^^ with ap G M. 

Theorem 16.31 provides an alternative to the classical method of rewriting a poly- 
nomial / in terms of orthogonal polynomials by substituting each monomial in / 
by applying Theorem 16.11 It gives linear rules to rewrite the elements of G and 
the remainder of a polynomial divided by G as linear combinations of orthogonal 
polynomials of low enough order. The proof is in Appendix [71 

Theorem 6.3. 

(1) Span(K) = Span(ClB); 

(2) Let G he the reduced a-Grdhner basis of I(T>). The polynomial g £ G with 
LT{g) = is uniquely written as 

g = '^a- ^ bpiTii 

where b — [6/3]^eL,^<„Q solves the linear system [^/siz)]^^^ ^^j^ ~ 
[7ra(z)]^gp; in words the coefficient matrix is the evaluation matrix over 
T> of the orthogonal polynomials with x^ in tail of g and the vector of 
constant terms is the evaluation vector of tTq . 

(3) Let p e M[x] be a polynomial and [p{z)]z£'d its evaluation vector. The 
polynomial p* defined as 

(14) P* 

where a = [ap]p^L solves the linear system [Trp{z)]^^^ ^^j^a = [p{z)]^^jj, 
is the unique polynomial such that p*{z) = p{z) for all d £ T> and p* G 
Span(C'6). 

Theorem 16.31 provides a pseudo-algorithm to compute a Grobner basis for I{T>) 
and interpolating polynomials at T> in terms of orthogonal polynomials of low order 
directly from T> and OB. In Table [3] we give a variation of the Buchberger-MoUer 
algorithm [7] which starting from a finite set of distinct points T> and a term- 
ordering a returns L and the expressions g = tTq — J^peL I3< a ^P^l^ ^o^' .9 f^ie 
reduced a-Grobner basis of I{T>). It does so by performing linear operations. If 
the real vector [p(z)]zgx' is assigned, then the expression p* = J^/si^l '^^"'0 '^^^ ^'^^ 
be found using Item 3 in Theorem 16.31 This permits to rewrite every polynomial 
p € W[x] as a linear combination of orthogonal polynomials. 

Summarising: given a function /, a finite set of distinct points T> C M.'^ and a 
term-ordering cr, a probability product measure A"* over M'*, its system of product 
orthogonal polynomials, and a random vector with probability distribution A'^, then 
the expected value of / with respect to A'' can be approximated by 

(1) computing L with the algorithm in Table [3] and 

(2) determining the unique polynomial p* such that p*{z) = f(z) for all z £ T>, 
by solving the linear system ['^piz)]^^jj a = [f{z)]j^^jy. The polynomial 
p* is expressed as linear combination of orthogonal polynomials. 

(3) The coefRcient oq of ttq is the wanted approximation. 

Recall that p*{x) — X^zep /(-^)^<i(^) is a linear combination of the indicator func- 
tions of the points in T> (Lagrange polynomials) and hence qq = X^zeX" •^(^) ^ i^d.{X j) 
In particular, ¥,^d {lz{X)) = X^, z Cz T>, can be computed by applying the above 
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Input: a set T) of distinct points in M'', a term-ordering a and any vector norm || • ||. 
Output: the reduced tr-Grobner basis G of I{'D) as linear combination of orthogonal 
polynomials and the set L. 
Step 1: Let L = {0 S Z^^), OB ^ [1], G = [ ] and M = [xi, . . . 
Step 2: If A/ = [ ] stop; else set a;" — mino-(M) and deleted from M . 
Step 3: Solve in h the overdetermined linear system [^^{z)]^^^ ^^^h = 
[7rQ,(z)]^gP and compute the residual 

P = - [^f}{z)]^^v,p&L ^ 

Step 4: 

(1) If IIpII > 0, then include a in L, and include in M those elements of 
{a;ia;", . . . , x^x"} which are not multiples of an element in M or of 
LT{g), g G G. Return to Step 2. 

(2) If IIpII = 0, then include in Gthe polynomial 

where the values bp, ^ d L, are the components of the solutions b of 
the linear system in Step 3. Delete from M all multiples of x". 
Table 3. Buchberger-MoUer algorithm using Orthogonal Polynomials 



to / = Iz- Notice however that as A'' is a product measure, the Xz can be ob- 
tained from the one-dimensional ones as noticed before Theorem 15.11 It would be 
interesting to generalise this section to non-product measures. 

6.1. Characterisation of polynomial functions with zero expectation. In 

this section we characterise the set of polynomials with the same expected value. 

As mentioned in Section 1 1.1 1 given V, its vanishing ideal 1(1^), a term-ordering 
a and a Grobner basis G of /(I?) with respect to a, then any polynomial p € R[x] 
can be written as 

geG 

where r{x) is unique in Span(S) such that r{z) =^ p{z) for all z G I? and can be 
written as r{x) = J2zev Pi^^dix), where Iz, z £ V, are the product Lagrange 
polynomials in Section [S] Theorem 16.31 savs how to write r over OB. 

If p G m.[x] is such that E {p{X)) = E {r(X)) then E {p{X) - r{X)) = 0. Fur- 
thermore we have p — r E I{T^). Hence instead of studying directly the set 

{pGK[a;] :E(p(X))=E(r(X))} 

we characterize the set 

8o^{gEl{V):K{g{X))^Q} 

in Theorem 16.41 whose proof can be found in Appendix [T) Hence if p G M\x\ is 
such that p = g + r with g E So and r G Span(S) then by linearity E(]3) = 

E(zi,...,.„)ei5M2i'---'^»)^"i 

Theorem 6.4. Let X"^ be a product probability measure with product orthogonal 
polynomials TTa{x), a G Z>q and let X be a random vector following A''. Let 
D C M.'^ be a set of distinct points, a a term- ordering, G be the a-reduced Grobner 
basis of I(T>) written as linear combination of orthogonal polynomials, that is for 
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g (z G and x"" = LT{g) write 

Let p = J^geG^a^ ^ /"'^ suitable Qg G M.[x], and consider the Fourier 

expansion of each Qg , g £ G, 

(15) Qg = C^fe)""/? +Co(gg)7I-ag + Cpi^a)'^P 

Then E {p{X)) — if and only if 

(16) J2\M\l^0iqg)^T. E Mlcpiqg)cp{g) = 
geG geGag>„l3eL 

Note that any polynomial can be written according to Equation (1151) where the 
first sum includes terms higher in a than LT{g) and terms that do not appear in 
g, the last sum includes terms lower than LT{g) in a. Hence, the key observation 
in the proof is that E(7rm 7r„) = is n 7^ m and then linearity of E is used. 
Importantly, only terms of low enough Fourier order (the second and third terms 
in Equation p5p l matter for the computation of the expectation. 

Example 6.5. Consider Zi and Z2 two independent standard normal random 
variables and hence the Hermite polynomials. Consider also the five point design 

V = {(-6, -1), (-5, 0), (~2, 1), (3, 2), (10, 3)} 

and the a =DegLex term-ordering over the monomials in R[a::,j/]. The algorithm 
in TableElgives OB = {1, Hiiy), Hi{x), Hiix)Hi{y), H2ix)} and G = {51,52,53} 
where 

gi^H2iy)-H,{x) + 2H,iy)-4 

52 = H2{x)Hi{y)-9H2{x)+47Hi{x)Hi{y) - 123Hi(a;) + 27mi(y) - 399 

53 = Hsix) - A7H2{x) + iQQHi{x)Hi{y) ~ 845i/i(a;) + 2040i7i(y) - 2987 

By Theorem 16.41 for the purpose of computing its expectation a polynomial p = 
9i5i + 9252 + 9353 G -^(2^) can be simplified to have the form 

p = (4^'i/2(y) + c[}|o)Hi(x) + ^ll^^HM + c[J^o))5i 
+ (4')i/2(:r)Hi(y) +cgo)^2(:r) + cjj|,)i/i(x)i/i(y) +c[^;„)i/i^ 
+ {cfH^{x) + c[^;o)^2(x) + S^l^^H^{x)H^{y) + S^l^^H^{x) + S^l^^HM + c^t^^Qz 

and furthermore by Equation p6p 

(1)21 +2c^'^ -4c^'^ 

'-O ^- '-(1.0) + ^'-(0,1) ^'-(o.o) 

I „(2)9| _ Q„(2) 91 I /I7^(^'' —19'?^^^'' ^971^'-^) _ '^QQ^'^^) 

+ c^^h\ - 47c[2jo)2! + 300c[J^) - 845c[iJ^o) + 2040cjJ^) - 2987c[o^o) ^ 

In practice, for i = 1, 2, 3, put coefficients of gi and qi in two vectors, multiply them 
component wise and sum the result. For example the above equation is satisfied by 



2 

5863/2987 





= -34 


cn 0) 


= 


^(1) 


= -2 


''(0,0) ^ 


-8 








(2) 

Co 


= 


(2) 
C(2.0) 


= 


^(2) 
''(1,1) 


= 1 


(2) 
C(l 0) - 


-2 


(2) 
C(0,1) 


= 1 


(2) 
C(0,0) ~ 


Cq 


= 10 


''(2,0) 


= 2 


'^(1,1) 


= 


c^-'^) - 
C(1,0) - 


-5 


C(0,1) 


= 1 


„(3) _ 
(0,0) ~ 
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by adding Hi{x) to qi and Hi{y) to (72 which are not influent in the computation 
of the expectation we get the foUowing zero mean polynomial 

p{Zi,Z2) = 

ZlZl + VdZl + Z\Zl - 9ZlZ^ + 47ZiZ| - 469Zf + 3002Z^Z2 + ZfZ^ - GZ^Z^ 

- 123^1 + 270Z| - 8614Z^ + 20990^1^^2 + 96^1^^^ - 282ZiZ| - 424Z| 

71 78^81 4 

- 87898560/2987^1^ - 6700ZfZ2 + 1389ZiZ| - 1690Z| + 

2987 

218275468^ ^ .„..„9 307862660 ^ 5937584 „ 5931425 



-Z1Z2 + 4845Z^ + """ Zi - """'""' Z 



2987 ' " ' ' 2987 ' 2987 ^ 2987 

7. Appendix: proof 

Theorem [gTH 

Proof. 1. The proof is by induction on the monomial degree fc. From the three 
terms recurrence formula tTj+i = (jjx — aj)iTj — jSjiTj^i we have 

XTTj — — \ -TT-i H -TTj-l 

I3 Ij I3 

For fc = we have cc" — t^o{x) = co(a;'^)7ro. For k — 1 from the three terms 
recurrence formula we have 

X — XTTo — \ VTo = ci{x)tti + co(a;)7ro 

70 70 

In the inductive step the thesis holds for k and we prove it for k + 1. From the 
three recurrence formula we have 

x'^' = xx' = V c,{x')x., = ^^■(-') + + ^^^-^ 



j=o j=o 



7j 7j 7i 



k fc-1 



i=i ^^"^ j=o 'J J=0 



fe-1 



h Cj(X ) h C-, + l(a; j TTj H TTfe H TTfc+l 

7j-i 7j 7j+i/ 7fc-i 7fc 



H TTO H TTfe H TTo 

70 7fe 71 

fc-i / Jk) (fc) \ / (fc) (fc) 



This concludes the proof of the first part of the theorem. To prove the second part 
we apply what we just proved and unfold the multiplication. 

Given = x"^ ■ ■ ■ x'^'' , the polynomial tTq, = tTq-j (xi) ■ ■ ■ tTq^ (xd) is the product 
of d univariate polynomials tTq,^ each of degree aj in Xj, j = 1,. . . ,d. Clearly if 
aj = then tt^j — 1 and Xj does not divide x". Furthermore we have 

-.=nE4s^ 

j = l k=0 

We deduce that tTq is a linear combination of x" and of the power products which 
divide x", that is of power products x^ with (3 < a component wise. Vice versa, 
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applying the first part of the theorem we have 



fc=i 



fe=i 



Jk=0 



and commuting product with sum shows that x" is a Unear combination of products 
of Trp.{xi) where /3 = (/3i, . . . , (3d) is such that /3 < a component wise, that is x^ 
divides x". □ 



Theorem [6T3l 

Proof. RecaU that B and OB are defined in terms of a common set L of d-dimensional 
vectors with non-negative integer entries satisfying the property of 'factor-closeness' 
, that is if (q!i, . . . , ad) G L then . . . , /3d) G L if and only if /3i < Ui for all 
i = 1, . . . , d. 

(1) If e S for some a, then by Theorem 16. II = X]/3<q ^^'''/J Since (3 < a 
then (3 ^ L and so each tt^ e OS: we have that cc" belongs to Span(C';B). 
The vice- versa is proved analogously. 

(2) The matrix [t: jj{z)\z<z-D ,p£L is a square matrix since L has as many elements 
as T) and has full rank. The linear independence of the columns of such 
a matrix follows from the fact that each linear combination of its columns 
corresponds to a polynomial in Span(C';S) which coincides with Span(S) 
whose elements do not vanish at T). 

Any polynomial g G G in the Gobner basis can be written as 

g = x°' - ^ Cfjx^ 

where x" = LT{g) is a multiple of an element of B. By Theorem 16.11 we 
have 

The polynomial tTq appears only in the first sum, for the other terms in the 
first sum observe that as /3 < a then (3 E L and also (3 <a- ct- Analogously, 
for the second sum we consider 7 < /3 < a; since (3 E L then 7 g i and 
since 7 < a then 7 a. And so, with obvious notation. 

Since g{z) = for z G 2?, then the vector b = [bp]p of the coeffi- 
cients in the identity above solves the linear system [TTf;{z)]z^x).a >cPeLb = 
[T^aiz)]zev- Furthermore, since [TTi3{z)]zev.ag>„i3€L is a full rank matrix, 
then 6 is the unique solution of such a system. 

(3) Let p* = X]/3g-L ^z?""/? be a polynomial whose coefficients are the elements 
of the solution of the linear system [Trf}{z)]zev,0eLO. ~ [p{z)]z^v- Such a 
polynomial obviously interpolates the values p{z), z E T), and, since the 
columns of [T^p{z)]zev,pe l are the evaluation vectors of the elements of OB 
at 2?, it belongs to Span(C'S). We conclude that p* is the unique polynomial 
interpolating the values p{d), d £ V, w.r.t. {V, Span(OS)). 

□ 

Theorem [6T4l 
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Proof. As G is a Grobner basis of /(I?), then for every p € I I'D) and g E G there 
exist Qg G R[x] such that p — J2geG 1g9- Since by hnearity 

\S6G / geG 

the thesis follows once we show that, for each g £ G and x"' = LT{g) 
^{(lg9) ^hajlcoiqg) - Y c^(99)c^<(5)lk/3liA 

holds. 

From Equation (|15p we have 

9s5= Y cp{qg)'Kpg + co{qg)TTo,^g+ ^ cp{qg)Trfsg 

and substitute the Fourier expansion of g given in Theorem 16.31 

In computing the expectation we use the fact that E {iThTTk) — for different mono- 
mials h and k. Then the expectation of the first sum vanishes, the expectation of the 
middle term gives UtTq^ \\lco{qg) and the last sum gives - I]ag>„/3eL <^l^ilg)ci3ig)\Ui3\\l 

□ 
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